Mixed effects models lme4

library(dplyr)
# For Model fitting
library(lme4)
library(purrr)
# For diagnostics
library(performance)
# For adding new columns
library(tibble)
# For highlighting min value in a table
library(reactablefmtr)
# For boxcox function
library(caret)

# Load data
sys.source("./codes/scripts/code_join_data_full_dataset.R", envir = knitr::knit_global())

# Load functions

## Models
sys.source("./codes/functions/functions_models.R", envir = knitr::knit_global())

## Inference
sys.source("./codes/functions/function_anova_and_tukey_tables.R", envir = knitr::knit_global())

## Plot
sys.source("./codes/functions/function_masomenos_plots.R", envir = knitr::knit_global())

Select performance and traits variables

data_for_models <- 
    data_complete %>%
        
        # Calculate nitrogen use efficiency column
        # I followed Leaf traits explaining the growth of tree 
        # species planted in a Central Amazonian disturbed area
         mutate(photo_nitrogen_use = (amax*sla_cm2_g*0.1)/perc_n) %>% 
        
        # select variables that are going to be used in the models
        select(id, spcode, treatment,nfixer, init_height, 
               
               # Plant preformance
               total_biomass, above_biomass, below_biomass, agr,
               
               # Traits
               amax, gs, wue, d13c, d15n, photo_nitrogen_use) %>% 
        mutate(id = factor(id))


data_for_models$nfixer <- relevel(data_for_models$nfixer,ref = "fixer" )
levels(data_for_models$nfixer)[1]

[1] "fixer"

Box-Cox transformation

# Nonnormality and heterogeneity of variance were corrected by applying  box-cox transformation
box_cox <- 
    data_for_models %>%
        
        # Transforme only the traits
        select(5,6:9,10:15) %>% 
        purrr::map(., ~ BoxCoxTrans(.x))

Add variables transformed to the dataset

data_for_models_transformed <- 
    cbind(data_for_models, 
      
        # Plant's performance
        total_biomass_transformed = predict(box_cox$total_biomass, data_for_models$total_biomass),
        above_biomass_transformed = predict(box_cox$above_biomass, data_for_models$above_biomass),
        below_bioomass_transformed = predict(box_cox$below_biomass, data_for_models$below_biomass),
        agr_transformed = predict(box_cox$agr, data_for_models$agr),
      
        # Traits
        amax_transformed = predict(box_cox$amax, data_for_models$amax),
        gs_transformed = predict(box_cox$gs, data_for_models$gs),
        wue_transformed = predict(box_cox$wue, data_for_models$wue),
        
        # d13 and d15 where not transformed because the data has negative values
        
        photo_nitrogen_use_transformed = predict(box_cox$photo_nitrogen_use, data_for_models$photo_nitrogen_use),
      
        # Covariate
        init_height = predict(box_cox$init_height, data_for_models$init_height))  %>% 
        
        # Remove original variables (non-transformed)
        select(id, spcode, treatment, nfixer, init_height, everything(),
               -c(5,6:9,10:12,15))    

Models: Questions 1 and 2

\[response\sim treatment*fixer\ + initial\ height + random( 1|\ specie)\]

# Take response variables' names 
response_vars_q1_q2 <- set_names(names(data_for_models_transformed)[6:(ncol(data_for_models_transformed))])

model_list_q1_q2 <- map(response_vars_q1_q2, ~ mixed_model_1(response = .x, data = data_for_models_transformed))

Models: Questions 3

\[performance\sim treatment*fixer*trait\ + initial\ height + random( 1|\ specie)\]

# This funtion takes two columns a create a model formula with all the possible 
# combinations

model_combinations_formulas <- function(y_var, x_var){
           
        variables <- crossing(y_var, x_var)
        pattern <- c('\\+.*|nfixer|treatment|[[:punct:]]| ')
        
           
        # Model
        models <- paste0(variables$y_var,"~treatment*nfixer*",
                            variables$x_var,"+init_height+(1|spcode)")
           
        formulas <- map(models, as.formula)
        
        # Add name to the list 
        names(formulas) <- stringr::str_replace_all(formulas, pattern, replacement = '')
    return(formulas)
}

Scale preditors

data_for_models_transformed_scaled <-         
    data_for_models_transformed %>%   
        
        # scale pedictors
        mutate(across(c(5, 10:15), scale))

# Select traits (x_vars)
traits_names <- set_names(names(data_for_models_transformed_scaled)
                          [c(6,7, 12:15)])

# Select plants performance vars (y_vars)
performance_names <- set_names(names(data_for_models_transformed_scaled)[8:11])


models_lmer_formulas <- model_combinations_formulas(performance_names, traits_names)

model_list_q3 <- map(models_lmer_formulas, 
                       ~ lmer(.x, data = data_for_models_transformed_scaled))

Validation plots performance package

Models for questions 1,2

map(model_list_q1_q2, check_model)

$d13c

$d15n

$total_biomass_transformed

$above_biomass_transformed

$below_bioomass_transformed

$agr_transformed

$amax_transformed

$gs_transformed

$wue_transformed

$photo_nitrogen_use_transformed

Models for question 3

map(model_list_q3, check_model)

$`abovebiomasstransformed~amaxtransformed`

$`abovebiomasstransformed~d13c`

$`abovebiomasstransformed~d15n`

$`abovebiomasstransformed~gstransformed`

$`abovebiomasstransformed~photonitrogenusetransformed`

$`abovebiomasstransformed~wuetransformed`

$`agrtransformed~amaxtransformed`

$`agrtransformed~d13c`

$`agrtransformed~d15n`

$`agrtransformed~gstransformed`

$`agrtransformed~photonitrogenusetransformed`

$`agrtransformed~wuetransformed`

$`belowbioomasstransformed~amaxtransformed`

$`belowbioomasstransformed~d13c`

$`belowbioomasstransformed~d15n`

$`belowbioomasstransformed~gstransformed`

$`belowbioomasstransformed~photonitrogenusetransformed`

$`belowbioomasstransformed~wuetransformed`

$`totalbiomasstransformed~amaxtransformed`

$`totalbiomasstransformed~d13c`

$`totalbiomasstransformed~d15n`

$`totalbiomasstransformed~gstransformed`

$`totalbiomasstransformed~photonitrogenusetransformed`

$`totalbiomasstransformed~wuetransformed`

Save lists with the models

saveRDS(model_list_q1_q2, file = "./processed_data/models_q1_q2.RData") 
saveRDS(model_list_q3, file = "./processed_data/models_q3_3_way_interaction.RData")